Modelos lineales

Regresión y factoriales (ANOVA y ANCOVA)

Organización del día

Hora

9:30 a 10:30 teoría
10:30 a 11:30 ejercicios
11:30 a 11:45 break
11:45 a 12:45 teoría
12:45 a 13:30 ejercicios y dudas

Temario

Modelos lineales regresión
Ejercicio regresión
Modelos lineales factoriales
Ejercicio ANCOVA

Introducción modelos lineales

Definición
Tipos
Supuestos
Flujo de trabajo

Definición

Relación lineal entre dos o más variables, dependiendo una (variable dependiente) de la(s) otra(s) (variable(s) explicativa(s))

Variable dependiente: y Variable(s) explicativa(s): x

Modelo matemático lineal: y = 2x o y = 3 - 2x

Objetivos

Descripción de la relación
Inferencia o generalización de los resultados
Predicciones

Tipos

Variable dependiente: siempre cuantitativa

Variables explicativas:

Cuantitativas: modelos de regresión (simple o múltiple)
Cualitativas o categóricas: ANOVA
Cuantitativas y categóricas: ANCOVA

Supuestos del modelo

Independencia
Linealidad
Normalidad
Homocedasticidad

Flujo de trabajo

Pasos a seguir análisis de datos con modelos lineales (Cayuela y de la Cruz, 2022)

Regresión simple y regresión múltiple

Estructura

Ecuación regresión:

Parte determinista: y_i = β₀ + β₁ x₁ + Ɛ_i

Parte estocástica: Ɛ_i ~ N(0,σ²)

β₀: ordenada en el origen

β₁: pendiente

Ɛ_i: residuos

Residuos

Ajuste modelo lineal

Función lm(y ~ x)

Ejemplo

Pregunta: Masa corporal especies pingüinos y cómo esta es influenciada por la longitud de la aleta de los pingüinos

Hipótesis

H₀ (hipótesis nula): no existe relación entre la masa corporal de los pingüinos y la longitud de sus aletas

H₀ : β₁ = 0

H₁ (hipótesis alternativa): pendiente de la recta que relaciona masa corporal con longitud aleta es distinta de cero

H₀ : β₁ ≠ 0

1. Leer los datos en R y explorar gráficamente la relación entre las variables

library(palmerpenguins)
data(penguins)
head(penguins)

# A tibble: 6 × 8
  species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
  <fct>   <fct>              <dbl>         <dbl>             <int>       <int>
1 Adelie  Torgersen           39.1          18.7               181        3750
2 Adelie  Torgersen           39.5          17.4               186        3800
3 Adelie  Torgersen           40.3          18                 195        3250
4 Adelie  Torgersen           NA            NA                  NA          NA
5 Adelie  Torgersen           36.7          19.3               193        3450
6 Adelie  Torgersen           39.3          20.6               190        3650
# ℹ 2 more variables: sex <fct>, year <int>

# eliminamos NAs variables que vamos a utilizar y pasamos masa corporal a kg
library(dplyr)
penguins_clean <- penguins %>%
  filter(!is.na(flipper_length_mm), !is.na(body_mass_g))
penguins_clean$body_mass_kg <- penguins_clean$body_mass_g/1000

library(ggplot2)
ggplot(data = penguins_clean, aes(x = flipper_length_mm, y = body_mass_kg)) +  
  geom_point(color = "blue", size = 2) +      
  labs(x = "Longitud aleta (mm)",         
       y = "Masa corporal (kg)",                   
      title = "Relación entre masa corporal y longitud aleta "   
  ) +
  theme_minimal()

2. Ajuste modelo de regresión

lm_reg1 <- lm(body_mass_kg ~ flipper_length_mm, data = penguins_clean)

3. Variación explicada por el modelo y resolución hipótesis

Tabla del análisis de la varianza

anova(lm_reg1)

Analysis of Variance Table

Response: body_mass_kg
                   Df  Sum Sq Mean Sq F value    Pr(>F)    
flipper_length_mm   1 166.453 166.453  1070.7 < 2.2e-16 ***
Residuals         340  52.855   0.155                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Suma de cuadrados: variabilidad de la variable y

Coeficiente de determinación (R²): bondad de ajuste del modelo (0-1)

4. Interpretación del modelo

summary(lm_reg1)


Call:
lm(formula = body_mass_kg ~ flipper_length_mm, data = penguins_clean)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.05880 -0.25927 -0.02688  0.24733  1.28869 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -5.780831   0.305815  -18.90   <2e-16 ***
flipper_length_mm  0.049686   0.001518   32.72   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3943 on 340 degrees of freedom
Multiple R-squared:  0.759, Adjusted R-squared:  0.7583 
F-statistic:  1071 on 1 and 340 DF,  p-value: < 2.2e-16

y = 0.05x - 5.78

Fiabilidad estimadores para predecir

R² y R² ajustado

5. Predicciones con el modelo

library(broom)
library(tidyverse)
library(dplyr)
preds <- augment(lm_reg1, type.predict = "response", se_fit = TRUE)
preds <- preds %>%
  mutate(lower = .fitted - 1.96 * .se.fit,upper = .fitted + 1.96 * .se.fit)
ggplot(data = preds, aes(x = flipper_length_mm, y = body_mass_kg)) +  
  geom_point(color = "blue", size = 2) +     
  geom_line(aes(y = .fitted)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2) +
  labs(x = "Longitud aleta (mm)",         
       y = "Masa corporal (kg)",                   
      title = "Predicciones modelo de regresión"   
  ) +
  theme_minimal()

6. Revisión supuestos del modelo

Supuestos de linealidad, normalidad y homocedasticidad

par(mfrow=c(2,2))
plot(lm_reg1)

Supuestos linealidad, normalidad, homocedasticidad

Supuesto normalidad

Supuesto homocedasticidad

Datos sobreinfluyentes en el análisis

Normalidad:

shapiro.test(residuals(lm_reg1))


    Shapiro-Wilk normality test

data:  residuals(lm_reg1)
W = 0.99301, p-value = 0.1123

Linealidad:

library(lmtest)
resettest(lm_reg1)


    RESET test

data:  lm_reg1
RESET = 15.998, df1 = 2, df2 = 338, p-value = 2.3e-07

Homocedasticidad:

library(car)
ncvTest(lm_reg1)

Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 2.177953, Df = 1, p = 0.14

Si no se cumplen supuestos del modelo: transformación variables (logarítmica, raíz cuadrada, Box-Cox, log-log, etc.)

Regresión múltiple

y_i = β₀ + β₁x₁ + … + β₁x_n + Ɛ_i

Ɛ_i ~N(0,σ²)

Suma de cuadrados tipo I o III

Suma de cuadrados: tipos (Cayuela y de la Cruz,) 2022

Suma de cuadrados tipo I: anova()

Suma de cuadrados tipo III: Anova(lm,type=“III”)(library(car))

Modelos factoriales: ANOVA y ANCOVA

Análisis de la varianza (ANOVA)

Hipótesis: comparación medias distintos grupos (niveles de factor)

H₀: medias de los grupos son iguales

H₀: µ_a = µ_b = … = µ

H₁: al menos una media es distinta a las demás

H₀: Ǝµ_j ≠ µ

Ejemplo

Pregunta: Diferencias ratio del pico (longitud/profundidad) de tres especies de pingüino

Hipótesis

H₀ (hipótesis nula): no existe diferencias en las medias del ratio del pico para las tres especies

H₀ : µ_adelie = µ_gentoo = µ_chinstrap = µ

H₁ (hipótesis alternativa): la media de al menos una de las tres especies es distinta a las demás

H₀ : Ǝµ_j ≠ µ

1. Explorar gráficamente la relación entre las variables

penguins2 = penguins %>%
  mutate(bill_ratio = bill_length_mm / bill_depth_mm)
ggplot(data = penguins2, aes(x = species, y = bill_ratio, fill = species)) +  
  geom_boxplot(alpha = 0.8, color = "gray30") +
  scale_fill_manual(values = c("#F8BBD0", "#C8E6C9", "#B3E5FC")) +
  labs(x = "Species",         
       y = "Ratio del pico",                   
      title = "Diferencias ratio del pico entre especies"   
  ) +
  theme_minimal()

2. Ajuste modelo factorial

penguins2 <- penguins2 %>%
  mutate(species = as.factor(species))
lm_reg2 <- lm(bill_ratio ~ species, data = penguins2)

3. Variación explicada por el modelo y resolución hipótesis

Tabla de análisis de la varianza

anova(lm_reg2)

Analysis of Variance Table

Response: bill_ratio
           Df Sum Sq Mean Sq F value    Pr(>F)    
species     2 75.766  37.883  1494.9 < 2.2e-16 ***
Residuals 339  8.591   0.025                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4. Interpretación del modelo

summary(lm_reg2)


Call:
lm(formula = bill_ratio ~ species, data = penguins2)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.60912 -0.10810  0.00319  0.10713  0.60467 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       2.11973    0.01295  163.62   <2e-16 ***
speciesChinstrap  0.53403    0.02325   22.97   <2e-16 ***
speciesGentoo     1.05587    0.01934   54.61   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1592 on 339 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.8982,    Adjusted R-squared:  0.8976 
F-statistic:  1495 on 2 and 339 DF,  p-value: < 2.2e-16

Variables indicadoras o variables dummy. Un nivel es la coordenada en el origen (nivel de referencia)

y_i,j = β₁x₁ + β₂x₂ + β₃x₃ + Ɛ_i

y_i = 2.12 + 0.53Chinstrap + 1.06Gentoo

Ver el nivel de referencia

 contrasts(penguins2$species)

          Chinstrap Gentoo
Adelie            0      0
Chinstrap         1      0
Gentoo            0      1

Comparaciones post-hoc entre los distintos grupos

library(multcomp)
posthoc_species <- glht(lm_reg2, linfct=mcp(species="Tukey"))
summary(posthoc_species)


     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lm(formula = bill_ratio ~ species, data = penguins2)

Linear Hypotheses:
                        Estimate Std. Error t value Pr(>|t|)    
Chinstrap - Adelie == 0  0.53403    0.02325   22.97   <2e-16 ***
Gentoo - Adelie == 0     1.05587    0.01934   54.61   <2e-16 ***
Gentoo - Chinstrap == 0  0.52184    0.02406   21.69   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

5. Representación gráfica diferencias significativas

ggplot(data = penguins2, aes(x = species, y = bill_ratio, fill = species, color = species)) +  
  geom_boxplot(alpha = 0.8, color = "black") +
  geom_jitter(width = 0.15, alpha = 0.6, size = 2, shape=21, color="black") +
  scale_fill_manual(values = c("#F8BBD0", "#C8E6C9", "#B3E5FC")) +
  scale_color_manual(values = c("#F8BBD0", "#C8E6C9", "#B3E5FC")) +
 labs(x = "Species",         
       y = "Ratio del pico",                   
      title = "Diferencias longitud aleta entre especies"   
  ) +
  theme_minimal() +
  geom_text(data = penguins2[2,], aes(y = 4, label = "a"), color="black", 
            size=5, nudge_x = 0) +
  geom_text(data = penguins2[2,], aes(y = 4, label = "b"), color="black", 
            size=5, nudge_x = 1) +
  geom_text(data = penguins2[2,], aes(y = 4, label = "c"), color="black", 
            size=5, nudge_x = 2)

6. Revisión supuestos del modelo

par(mfrow=c(2,2))
plot(lm_reg2)

Cambiar nivel de referencia del factor

levels(penguins2$species)

[1] "Adelie"    "Chinstrap" "Gentoo"

penguins2$species <- relevel(penguins2$species, ref = "Gentoo")
levels(penguins2$species)

[1] "Gentoo"    "Adelie"    "Chinstrap"

Análisis de la covarianza (ANCOVA)

Hipótesis: homogeneidad de pendientes

H₀: las pendientes son iguales

H₁: al menos una pendiente es distinta a las demás

Ejemplo

Pregunta: La masa corporal de los pingüinos se ve influenciada por la longitud de su pico, teniendo en cuenta el sexo del individuo

Hipótesis

H₀: no hay relación entre la masa corporal y la longitud del pico

H₀: descontando el efecto de la longitud del pico (si existe), no hay relación entre la masa corporal y el sexo del pingüino

Homogeneidad de pendientes

Asunción adicional a las de ANOVA y regresión: relación entre la variable respuesta y la covariable (pendiente de la regresión) es la misma en todos los grupos (homogeneidad de pendientes)

ggplot(subset(penguins, !is.na(sex)), aes(x = bill_length_mm, y = body_mass_g)) +
  geom_point(color = "steelblue", size = 2) +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  facet_wrap(~ sex) +
  labs(
    title = "Relación longitud pico y masa corporal",
    x = "Longitud pico (mm)",
    y = "Masa corporal (g)"
  ) +
  theme_minimal()

Ajuste modelo

lm_3 <- lm(body_mass_g ~ bill_length_mm*sex, data = penguins)

Comprobación de hipótesis

anova(lm_3)

Analysis of Variance Table

Response: body_mass_g
                    Df    Sum Sq  Mean Sq  F value    Pr(>F)    
bill_length_mm       1  74792533 74792533 191.8812 < 2.2e-16 ***
sex                  1  12051930 12051930  30.9194 5.572e-08 ***
bill_length_mm:sex   1    175714   175714   0.4508    0.5024    
Residuals          329 128239489   389786                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lm_4 <- lm(body_mass_g ~ bill_length_mm+sex, data = penguins)
anova(lm_4)

Analysis of Variance Table

Response: body_mass_g
                Df    Sum Sq  Mean Sq F value    Pr(>F)    
bill_length_mm   1  74792533 74792533 192.201 < 2.2e-16 ***
sex              1  12051930 12051930  30.971 5.427e-08 ***
Residuals      330 128415203   389137                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ggplot(subset(penguins, !is.na(sex)), aes(x = bill_length_mm, y = body_mass_g, color = sex)) +
  geom_point(size = 2) +
  scale_color_manual(values = c("#F8BBD0", "#B3E5FC")) +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(
    title = "Relación longitud pico y masa corporal",
    x = "Longitud pico (mm)",
    y = "Masa corporal (g)",
    color = "Sexo"
  ) +
  theme_minimal()

ggplot(data = subset(penguins, !is.na(sex)), aes(x = sex, y = body_mass_g, fill = sex, color = sex)) +  
  geom_boxplot(alpha = 0.8, color = "black") +
  geom_jitter(width = 0.15, alpha = 0.6, size = 2, shape=21, color="black") +
  scale_fill_manual(values = c("#F8BBD0", "#B3E5FC")) +
  scale_color_manual(values = c("#F8BBD0", "#B3E5FC")) +
 labs(x = "Sexo",         
       y = "Masa corporal (g)",                   
      title = "Diferencias masa corporal entre sexo"   
  ) +
  theme_minimal() +
  geom_text(data = penguins[2,], aes(y = 6500, label = "a"), color="black", 
            size=5, nudge_x = 0) +
  geom_text(data = penguins[2,], aes(y = 6500, label = "b"), color="black", 
            size=5, nudge_x = 1)

Pendientes significativamente distintas

Relación masa corporal con longitud de la aleta, teniendo en cuenta la especie

ggplot(subset(penguins, !is.na(species)), aes(x = flipper_length_mm, y = body_mass_g)) +
  geom_point(color = "steelblue", size = 2) +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  facet_wrap(~ species) +
  labs(
    title = "Relación longitud aleta y masa corporal",
    x = "Longitud aleta (mm)",
    y = "Masa corporal (kg)"
  ) +
  theme_minimal()

datos_limpios <- penguins %>%
  filter(!is.na(flipper_length_mm))
lm_5 <- lm(body_mass_g ~ flipper_length_mm*species, data = datos_limpios)
anova(lm_5)

Analysis of Variance Table

Response: body_mass_g
                           Df    Sum Sq   Mean Sq  F value    Pr(>F)    
flipper_length_mm           1 166452902 166452902 1211.946 < 2.2e-16 ***
species                     2   5187807   2593904   18.886 1.686e-08 ***
flipper_length_mm:species   2   1519564    759782    5.532  0.004327 ** 
Residuals                 336  46147424    137344                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(lm_5)


Call:
lm(formula = body_mass_g ~ flipper_length_mm * species, data = datos_limpios)

Residuals:
    Min      1Q  Median      3Q     Max 
-911.18 -251.93  -31.77  197.82 1144.81 

Coefficients:
                                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)                        -2535.837    879.468  -2.883  0.00419 ** 
flipper_length_mm                     32.832      4.627   7.095 7.69e-12 ***
speciesChinstrap                    -501.359   1523.459  -0.329  0.74229    
speciesGentoo                      -4251.444   1427.332  -2.979  0.00311 ** 
flipper_length_mm:speciesChinstrap     1.742      7.856   0.222  0.82467    
flipper_length_mm:speciesGentoo       21.791      6.941   3.139  0.00184 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 370.6 on 336 degrees of freedom
Multiple R-squared:  0.7896,    Adjusted R-squared:  0.7864 
F-statistic: 252.2 on 5 and 336 DF,  p-value: < 2.2e-16

y_i = -6787.3 + 4251.4 Adelie + 3750.1 Chinstrap + (54.6 - 21.8 Adelie -20.0) flipper_length

Gráfica de predicciones

newdata <- expand.grid(
  flipper_length_mm = seq(min(datos_limpios$flipper_length_mm), max(datos_limpios$flipper_length_mm), length.out = 100),
  species = levels(datos_limpios$species)
)
predicciones_grid <- augment(lm_5, newdata = newdata)

ggplot(datos_limpios, aes(x = flipper_length_mm, y = body_mass_g, color = species)) +
  geom_point(alpha = 0.5) +
  geom_line(data = predicciones_grid, aes(x= flipper_length_mm, y = .fitted, group = species, color = species), size = 1.2) +
  scale_color_manual(values = c("#F8BBD0", "#C8E6C9", "#B3E5FC")) +
  labs(
    title = "ANCOVA: predicciones por especie",
    x = "Longitud aleta (mm)",
    y = "Masa corporal (g)",
    color = "Especie"
  ) +
  theme_minimal()

Supuestos del modelo

par(mfrow=c(2,2))
plot(lm_5)